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Abstract 

Despite the success of the combination of conservative schemes and stag¬ 
gered constrained transport algorithms in the last hfteen years, the accurate 
description of highly magnetized, relativistic flows with strong shocks repre¬ 
sents still a challenge in numerical RMHD. The present paper focusses in the 
accuracy and robustness of several correction algorithms for the conserved 
variables, which has become a crucial ingredient in the numerical simula¬ 
tion of problems where the magnetic pressure dominates over the thermal 
pressure by more than two orders of magnitude. 

Two versions of non-relativistic and fully relativistic corrections have been 
tested and compared using a magnetized cylindrical explosion with high mag¬ 
netization (^10^) as test. In the non-relativistic corrections, the total energy 
is corrected for the difference in the classical magnetic energy term between 
the average of the staggered helds and the conservative ones, before (CAl) 
and after (CAl’) recovering the primitive variables. These corrections are 
unable to pass the test at any numerical resolution. The two relativistic ap¬ 
proaches (CA2 and CA2’), correcting also the magnetic terms depending on 
the flow speed in both the momentum and the total energy, reveal as much 
more robust. These algorithms pass the test succesfully and with very small 
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deviations of the energy conservation and very low valnes of the 

total momentnm In particnlar, the algorithm CA2’ (that corrects 

the conserved variables after recovering the primitive variables) passes the 
test at all resolntions. 

The nnmerical code nsed to rnn all the test cases is briefly described. 

Keywords: nnmerical RMHD, conservative schemes, constrained-transport 
schemes 


1. Introduction 


The popularization of the n se of conserya.tive methods in numerical RMHD , 
intro duced in the late nineties flKoide et ahl.ll999l:lKomissarovl.ll999al:lBalsaral. 
20011), has revolutionized the research in several helds of Relati vistic Astro¬ 
physics as, e.g., the study of the jets emanating from AGN (IKomissarovl. 


2003 


1999bl: ILeismann et ahl. l2005l: iKeppens et ahl. 120081: iMignone at ah 


2012 ) ^ the structure and dynamics of pulsa r wind nebulae (iKomissarov fc LvubarskvI . 


2010 ; 


Miznno et ah 


Del Zanna et al.l. l2004i IPorth et ahl. l2013h. or the production of GRB 


I Komissarov et al. . 20091 2010 : Rezzolla et al. . 2011 ). 


To prevent the generation of artihcial forces that can falsify the so¬ 
lutions, conservative schemes in numerical RMHD must be supplemented 
with some additional procedure to keep the magnetic held solenoidal along 
the simulation. Among the most succesful strategies is the combination of 
numerical algorithms in conservation form to evolve cell-centered (hnite- 
difference, hnite-volume) representations of the hydrodynamical variables, 
and constrained transport (CT) algorithms for the magnetic held discretized 
on a staggered grid. The CT scheme, based in the conservation of mag¬ 
netic hux across closed surfaces, maintains the value of the divergence of 
the magnetic held to the accuracy of machine round-oh errors on a spe- 
cihc discretization. It was o rig inally developed for artih c ial viscosity meth - 


ods ( Evans fc Hawlevl . 


1988 ). iDai fc WoodwardI (119981) . iRvu et al.l (119981) . 


and iBalsara fc Spicerl (119991) combined the CT disc retization with conser¬ 


vative schemes. Londrillo fc Del Zanna ( 2000 . 20041) developed the upwind 


constrained transport (UCT) strategy, which extends the C T method to 


1 999al 

Del Zanna et ah. 

2005: 

Anton et ah. 

2006; 

2007; 

Del Zanna et ah. 2 


A t present, many RMHD codes (IKomissarov . 


2003: Leismann et ah 


Mignone fc Bodol.l2006l: 


20051: IShibata fc Sekiguchil. 


Mignone et al.l. 12007: ICiacomazzo fc Rezzolla 
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use this approach. 

The combination of conservative algorithms and staggered CT schemes 
has proven to be very succesful in the simulation of highly magnetized, rel¬ 
ativistic flows with strong shocks. However the existence of two sets of vari¬ 
ables dehned on different grids requires in the most extreme cases some ad¬ 
ditional work to m ake the conservative and CT steps ful l y consistent . The 


proposed solutions 

IKomissarovl. 1999a: 

Balsara & Soicer. 

1999: 

Toth. 

2000: 

Mimone & Bodo. / 

i006tl relv on correcting the conserved quantities after each 


time step to make them consistent with the staggered helds. The present pa¬ 
per explores the performance of several correction algorithms focussing in 
their accuracy and robustness. 


The evolution system guarantees the fulhllment of this constraint for an 
initially divergence-free magnetic held at all later times, but to satisfy the 
constraint in numerical simulations of MHD hows poses a challenge. 


2. Conservative methods for RMHD and Constrained Transport 

2.1. The RMHD equations in conservation form 

The equations of ideal RMHD, representing the conservation of rest-mass, 
momentum and energy of the magnetized huid together with the induction 
equation can be written as a system of conservation law^. In Minkowski 
spacetime and Cartesian coordinates {{i,j, k} = {x, y, z}) this system read^ 


dU dF*(U) 
dt dx^ 


( 1 ) 


where the vector of conserved variables, U, and the huxes, F*, are the column 
vectors. 


/ D \ 
E 

\b^ J 


( 2 ) 


'^The solenoidal condition for the magnetic field is guaranteed by the evolution system 
however it is by no means trivial to satisfy numerically. 

^Throughout this paper, besides using units in which the speed of light is set to unity, 
a factor \/47r is absorbed in the definition of the magnetic field. 
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( 3 ) 


/ Dv^ \ 

j _ — V /W 

Ev^ + p*v^ — b^B'^/W 

^ _ ykj^t ^ 

In these equations, D, S\ E and B^ are the rest-mass density, the momentum 
density of the magnetized fluid in j-direction, the total energy density, and 
the magnetic held measured in the laboratory frame, 

D = pW, (4) 

= ph*W^v^ ( 5 ) 

E = ph*W^ -p* - {by, (6) 

where p is the proper rest-mass density of the huid, p* is the total pressure, 
and h* is the specihc enthalpy including the contribution of the magnetic 
held. These two last quantities are dehned according to 

b^ 

P* =P+- (7) 

h* = 1 + e + p/p+ b‘^/p, (8) 

where p is the huid pressure and e its specihc internal energy, b^ {p = 
0,1, 2, 3) are the components of the 4-vector representing the magnetic held 
in the huid rest frame and 6^ stands for 6^6^, where summation over repeated 
indices is assumed, w* are the components of the huid 3-velocity in the 
laboratory frame, which are related to the how Lorentz factor, W, according 
to: 


W 


1 

\/l — V^Vi 


(9) 


The following relations hold between the components of the magnetic 
held 4-vector in the comoving frame and the three vector components B^ 
measured in the laboratory frame: 


6° 

6* 


ITB ■ V , 

TDi 


4 


( 10 ) 

( 11 ) 



where v and B denote the 3-vectors (n®, v^) and {B^, B^, B^), respectively. 

The square of the modulus of the magnetic field can be written as 




B^ 


+ (B ■ v)2 


( 12 ) 


with B^ = B^Bi. 

An equation of state that relates the thermodynamic variables, e.g., p = 
p{p,e), is needed to close the system. 

To make more clear the corrections to be performed on them, the momen¬ 
tum and total energy densities are written making explicit their dependence 
on the magnetic field in the laboratory frame, B 


= [phW^ + B^y - (v ■ B)B\ 

(13) 

2^2 _ / . b )2 

(14) 

E = phW^-p+ — + -^ 




2.2. Conservative methods for RMHD and Constrained Transport 

Conservative methods exploit the conservation properties of the system of 
equations and can be directly applied to solve the equations of ideal RMHD 
written in conservation form, Eq. O- 
In these methods. 


where is the value of the conserved variable U at the point {xi,yj,Zk) 

(hnite-difference approach) or its volume average at the cell centered in that 
point (hnite-volume approach), and 



C 


X 


1 

Axi 


■ptc _ -px 


(16) 


(and similar expressions for Bfj if). Quantities 'Fi_^.if 2 jk 

merical fluxes. In the Godunov-type methods (hnite-volume approach), these 
huxes are obtained from the solution of Riemann problems at cell interfaces 
{xi+i/ 2 ,yj, Zk), where the initial left and right states, \JL,i+i/ 2 ) Ui{,j+i /2 are 
reconstructed values from the corresponding cell averages. 
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The algorithm dehned in Eq. fflB]) can be used to evolve system ([T]) in 
time. However, additional care has to be taken to prevent the divergence 
of the magnetic held to grow with time. Staggered CT algorithms main¬ 
tain the divergence of B exactly (i.e., to the accuracy of machine round-off 
errors) in each numerical cell. In this approach, the normal component of 
the surface-centered magnetic held on the cell interface at (a:i+i/ 2 ,2/j, is 
evolved according to a discretized version of the induction equation 


Ay ^^^+V2j+l/2,fc ^i+l/2j-l/2,fc) 

^(^i+l/2j,fc+l/2 “ ^i+l/2,j,fc-l/2)’ 

where quantities fl* are discretized representations of the components of f2 = 
vxB computed at cell edges in terms of spatial and temporal interpolations of 
the magnetic held and the velocity, or the numerical huxes of the conservative 
step. 

Once the staggered magnetic helds have been computed, the correspond¬ 
ing cell-centered helds can be obtained by interpolation. For second order 
accuracy, a linear interpolation is enough and 


'^^i+l/2,j,k 

dt 


B 


X 


2 i^i-l/2,j,k 


+ Bi+l/2,j,k) 


(18) 


(and similar expressions for Bfjj^ and B^j f^). 

In its simplest version, the evolution scheme consists of a conservative 
step in which the rest-mass, momentum and energy densities are advanced 
in time, and a CT step to advance the staggered magnetic helds and obtain 
the cell-centered ones. Finally, inherent to all the conservative methods in 
RMHD is to solve an implicit algebraic system to recover the primitive vari¬ 
ables (needed to compute the numerical huxes), V = {p,p,v\ B^) from the 
conserved ones, U, at every time step. 


2.3. Setting the problem 

According to the previous algorithm, at time t = the conserved 

variables have been advanced consistently with the magnetic helds dehned at 
cell centers in the previous time step and are hence consistent with a magnetic 
held at cell centers, ^cdu]k computed with the cell-centered 

scheme (conservative step). However these helds are diherent from those 
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defined at cell centers, ^stag!ij,fc> as the average of the staggered 

helds using Eq. ([TK]) . The fact that these two sets of magnetic helds are 
different is on the basis of the inconsistency of the algorithm, which makes 
it unsuitable for problems where the magnetic pressure dominates over the 
thermal pressure by more t han two orders of magnitude. _ 


The proposed solutions 

<omissarov. 1999al: Balsara & Spicer. 

1999; 

Toth 

2000; 

Mienone & Bodo. 

200 

3) rely on redehning the conserved quantities af- 


ter each time step to make them consistent with the staggered held. 


3. Correction of the conserved variables 


Rel ying in the approach pro posed bv lBalsara fc Spicerl fjlQQQIl for classical 
MHD, Mignone fc Bodo f 2006[l proposed a non-relativistic correction of the 
conserved energy after each time step according to 


h^stag .^cell 




d2 

■^stag 


(19) 


In the previous expression, Bceu and Bstag are, respectively, the magnetic 
helds at some cell as computed in the conservative step, and by averaging 
the staggered helds. E^ew and Estag are the corresponding conserved energy 
densities. 

The procedure is computationally efficient since it does not involve extra 
calls to the primitive recovery procedure per time step, although as the rest 
of correction algorithms proposed, it forces to evolve the two equations of the 
magnetic held components normal to the sweep in the conservative stepEI. On 
the other hand, the non-relativistic nature of the correction dismisses all the 
magnetic terms depending on the how velocity in the momentum and energy 
densities, making this correction questionable in those cases combining highly 
relativistic how velocities and large magnetic helds. 

This correction can be applied in two ways depending whether the correc¬ 
tion of the energy is done before (correction algorithm CAl) or after (CAl’) 


^ Since all the procedures are based on the advance in time of the cell-centered magnetic 
fields with the conservative algorithm, the number of three-dimensional variables increases 
in three. However, besides the additional calls to the primitive recovery procedure, where 
appropriate, the only extra computational cost is the pure advance of the two magnetic 
field components normal to the direction of the sweep, since these components are also 
reconstructed and used in the solution of the Riemann problems even in the case of no 
correction. 
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recovering the primitive variables. A potential drawback of the second ap¬ 
proach is that hydrodynamic, primitive variables are not made consistent 
with the staggered helds. 

In this work, we propose a full-relativistic correction in which the non- 
relativistic procedure CAl is used as a hrst approximation. The correction 
(CA2) proceeds as follows: 


1. Obtain an approximation to the total energy consistent with the stag¬ 
gered helds: 


r-'(l) _ rp 

■^stag .^cell 


r)2 _ o2 

-^cell -^stag 


( 20 ) 


2. Obtain an approximation to the primitive variables, namely from 
{L>, S, Bstag}- 

3. Use the how velocity v^) to complete the relativistic correction of the 
momentum and energy densities: 


s:.2 = s‘- (Bill - Bl^)v 


( 21 ) 


^(2) _ p(l) _ 

■^stag -^stag 


(t.0))2 


, o2 o2 1 , (''“’ ■ B»n)= - (v'O . B...j)= 

l-°cell -°stagJ T ^ 


4. Obtain the primitive variables from {D, Sg^ag; -^stag) Bstag}- 


( 22 ) 


Steps 1. and 2. correspond to the classical correction CAl. Steps 3. and 
4. complete the relativistic correction. The procedure involves one extra call 
to the recovery of primitive variables. 

Finally, as in the correction algorithm CAl’, the recovery of the primitive 
variables can be done prior to the correction of the conserved variables, lead¬ 
ing to a very simple and computationally efficient algorithm (CA2’). In this 
case, the primitive variables are recovered using the cell-centered magnetic 
helds as advanced by the conservative algorithm, Bceii, and then the con¬ 
served variables recalculate d from these primit ive variables and Bgtag- This 
is the approach followed by Komissarov f l999ah . As in the case of correction 
CAl’, a potential drawback of this approach is that hydrodynamic, primitive 
variables are not made consistent with the staggered helds (but, in contrast 
to CAl’, the corrected conserved variables do include the relativistic correc¬ 
tions). 












4. Numerical tests 


4 . 1 . The numerical code 

The basic ingredientes of the code used to test the correction procedure 
are the following: 


i) Cell reconstruction: second-order accurate values of the primitive vari¬ 
ables V at the left and right ends of the cells are obtained with linear 
functions and several limiters (MINMOD, VAN LEER, MC). MC and 
VAN LEER limiters can be degraded to MINMOD in case of strong 
shocks. No jump is allowed in the normal component of B at a cell 
boundary and the corresponding staggered magnetic field is used. 


ii) Riemann solvers: intercell numerical fluxe s are computed by means of 
HLL and HLLC ( Mignone fc Bodo . 20061) Riemann solvers. Accurate 
bounds of the maximum speeds of left and right propagating waves are 
obtained by solving the corresponding characteristic equation for the 
left and right states of each numerical interface. 


Ill 


V 


Time advance: the multidimensional equations of RMHD are advanced 
in time in an unsplit manner u sing TVD-preserving Run ge-Kutta meth¬ 


ods of second and third order fjShn fc Osherl . ll988l.ll989l ). The time step 

CFL 


is determined according to At = 




X mm 
hj 


/ Ax 

Ay 


\^y,hj\ 


(2D, pla¬ 


nar symmetry version), where Xx,ij and Xy^ij are the speeds of the 
fastest waves propagating in cell i,j along the x and y direction, re¬ 
spectively. 


iv) Constrained transport scheme as in lBalsara fc Spicer! fjl999l) 


Primi tive variables are recovered as in the TDy/ method of iNoble et ah 
(120061) and solving the resulting equation in Z = phW"^ by bisection. 


The code advances the total energy density without the rest-mass energy 
density, i.e., E — D. This strategy improves the performance of the conser¬ 
vative scheme when the total energy is dominated by the rest-mass energy. 
However, in this case it does not produce any effect, since in the selected 
test (see next Section) the total energy density is larger than the rest-mass 
energy density by several orders of magnitude. 
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4-2. The test: cylindrical magnetized bla st wave 

The setup for this test is taken from iKomissarovi (jl999al) . A cylindrical 
region of radius r = 0.8 with density p = 10“^ and thermal pressure p = 1 
is embedded in a static uniform medium with p = 10“^ and p = 3 x 10“®. 
A linear smoothing function is applied for 0.8 < r < 1. The whole region is 
threaded by a constant horizontal magnetic held in the x-direction, = 1. 
An ideal gas equation of state with 7 = 4/3 is used. 

The difference in pressure between the cylindrical region and the ambi¬ 
ent medium produces the expansion of the central region delimited by a fast 
forward shock propagating radially at almost the speed of light. Because 
of the strong sideways magnetic conhnement an elongated structure devel¬ 
ops in the x direction with a maximum Lorentz factor of kb ~ 4.5. This 
problem is particularly challenging because of the very large magnetization 
13 = 6V(2p) = 1.67 X lOT 

The test is solved in Cartesian (x, y) coordinates. The numerical grid 
covers a [— 6 , 6 ] x [— 6 , 6 ] square with the center of the cylindrical region at 
(0,0). Open boundary conditions are placed along the boundaries of the 
computational domain. 

A reference model using the HLLC Riemann solver and the third order 
Runge-Kutta (RK3) for time advance has been ran. For the cell reconstrucion 
we have used the VAN LEER limiter degraded to MINMOD (VLMM) when 
the relative jump in thermal or magnetic pressure within a shock exceeds 
0.5. The relativistic correction proposed in this paper (CA2) was chosen to 
correct the conserved variables after each time substep. A CFL of 0.45 and 
a numerical resolution of x Ny = 512^ cells complete the initial setup. 
Figured] shows the distributions of proper rest-mass density, gas pressure and 
magnetic pressure (in logarithmic scale), and flow Lorentz factor at t = 4.0 
as computed in the reference run. The fast magnetosonic shock and the 
elongated horizontal structure are clearly seen. 

Table [T] displays the total momentum and the relative change in total 
energy at t = 4.0 for the different correction algorithms and several numerical 
resolutions from 128^ to 1024^ cells. For the reference model (CA2, N^xNy = 
512^), the relative change in total energy is 2.60 x 10“®, whereas the total 
momentum (initially equal to zero) reaches a value of 4.11 x 10“^^. The 
same run without energy nor momentum corrections (CAO) crashes at the 
sixth iteration {t ~ 0.03), whereas the run u sing the cla. s sical c orrection 
(CAl) crashes at f ~ 0.39. The same run using |KomissaxQy| fjl999al) ’s recipe 
(CA2’), completes the test normally and with a relative change in total energy 
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Log(p) 



Figure 1: Proper rest-mass density, gas pressure and magnetic pressure (in logarithmic 
scale), and flow Lorentz factor at t = 4.0 as computed in the reference run of the cylindrical 
magnetized blast wave test discussed in the text. The magnetic field lines are printed on 
top of the Lorentz factor plot, logp G [-5.30,-2.71], logp G [-4.65,-0.94], \ogpm G 
[-0.60,-0.20], W G [1.00,4.63]. 
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Table 1: Total momentum (first row in two-row table entries) and relative change in total 
energy (secod row) in the cylindrical magnetized blast wave test discussed in the text at 
t = 4.0, as a function of the correction algorithm and numerical resolution. Table entries 
with a single number display the time at which the run stopped abnormally. 


A'x X Ny 

CAO 

CAl 

CAl’ 

CA2 

CA2’ 

128^ 

2.66 X 10-2 

1.13 

1.18 

1.19 X 10-^1 
1.51 X 10-^ 

8.00 X 10-12 
1.55 X 10-^ 

256^ 

2.29 X 10-2 

0.70 

0.73 

2.17 X 10-^1 
7.69 X 10-5 

7.50 X 10-1° 
7.91 X 10-5 

512^ 

2.72 X 10-2 

0.39 

0.42 

4.11 X 10-^1 
2.60 X 10-5 

1.55 X 10-® 
2.70 X 10-5 

10242 

3.75 X 10-2 

0.31 

0.37 

3.43 

7.94 X 10-° 
8.08 X 10-5 


of 2.70 X 10 ® and a total momentum of 1.55 x 10 The parameters chosen 
(VLMM,HLLC,RK3) produce a version of the code with very low dissipation, 
which makes the algorithms without correction algorithm (CAO) and with 
the classical correction (CAl) to fail in passing the test at any resolution. 
Code versions using corrections CA2 and CA2’ pass the test succesfully and 
with very small deviations of the energy conservation (;^10“^), and very 
low values of the total momentum (^10“®). In the case of the relativistic 
correction CA2, the run stopped abnormally for the largest resolution (1024^) 
close to the end of the test, at f = 3.43. 

In flows with high magnetization, accuracy problems can lead to unphys¬ 
ical internal energies and pressures. The correction algorithm CA2’ can be 
affected by this failure since in this case the recovery of primitive variables is 
done before the correction of the conserved variables and are then obtained 
from cell-centered representations of the magnetic held as advanced in the 
conservative step. In Table [21 the maximum and average relative differences 
in thermal pressure at f = 4.0 for the cases using the correction algorithms 
CA2 and CA2’ and numerical resolutions x Ny = 128^,256^,512^, are 
shown. The values of the average relative pressure difference between both 
algorithms tend to zero with numerical resolution, whereas the maximum rel- 
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Table 2: Maximum and average relative differences in thermal pressure for the cases 
using correction algorithm CA2 and CA2’, in the cylindrical magnetized blast wave 
test discussed in the text at t = 4.0. The maximum relative difference is computed 


as maxij 


bcA2(*, j) -PCA2'(*, j)| 


PCA2{i,j) 


The average relative difference is computed as 


Ezj bcA2(i,j) -PCA2'(b j)l 


Ei.,PCA2(bj) 


. In these expressions, i,j span the whole numerical grid. 


X Ny 

max 

aver 

128^ 

1.32 X 10-2 

5.53 X 10-® 

256^ 

2.77 X 10-2 

2.75 X 10-® 

512^ 

3.65 X 10-2 

8.65 X 10-® 


ative difference seems to increase with numerical resolution although keeping 
small values (si 1 — 3%). The conclusion is that the pressure values obtained 
with the algorithm CA2’ are consistent with those obtained with the algo¬ 
rithm CA2, in which the primitive variables are recovered once the conserved 
variables have been corrected. 

The same test with = 10,100 (magnetizations 1.67 x 10®, 1.67 x 10®, 
respectively) has been run without problems with the same algorithm choice 
(VLMM,HLLC,RK3) and correction algorithms CA2 and CA2’ with a spatial 
resolution of 512^ cells (see Fig. |2]). 

5. Summary and conclusions 

Despite the success of the combination of conservative schemes and stag¬ 
gered constrained transport algorithms in the last fifteen years, the accurate 
description of highly magnetized, relativistic flows with strong shocks flows 
represents still a challenge in numerical RMHD. The present paper focusses 
in the accuracy and robustness of several correction algorithms for the con¬ 
served variables, which has become a crucial ingredient in the numerical sim¬ 
ulation of problems where the magnetic pressure dominates over the thermal 
pressure by more than two orders of magnitude. 

Two versions of non-relativistic and fully relativistic corrections have been 
tested and compared using a magnetized cylindrical explosion with high mag- 
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Log(rho) 


Log(p) 



Figure 2: Proper rest-mass density, gas pressure and magnetic pressure (in logarithmic 
scale), and flow Lorentz factor at t = 4.0 as computed in the cylindrical magnetized blast 
wave test with = 100 (magnetization 1.67 x 10®) using the same algorithm choice as in 
the test run (VLMM,HLLC,RK3), correction algorithm CA2 and a spatial resolution of 
512^ cells. The magnetic field lines are printed on top of the Lorentz factor plot, logp € 
[-4.83,-2.75], logp e [-4.52,-1.00], logp^ S [3.698950,3.698986], IT e [1.00,4.68]. 
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netization (^10^) as test. In the non-relativistic corrections, the total energy 
is corrected for the difference in the classical magnetic energy term between 
the average of the staggered helds and the conservative ones, before (CAl) 
and after (CAl’) recovering the primitive variables. These corrections are 
unable to pass the test at any numerical resolution. The two relativistic ap¬ 
proaches (CA2 and CA2’), correcting also the magnetic terms depending on 
the flow speed in both the momentum and the total energy, reveal as much 
more robust. These algorithms pass the test succesfully and with very small 
deviations of the energy conservation and very low values of the 

total momentum (^10“®). In particular, the algorithm CA2’ (that corrects 
the conserved variables after recovering the primitive variables) passes the 
test at all resolutions. 

The numerical code used to run all the test cases is briefly described. 
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